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Three-dimensional stability of magnetically confined 
mountains on accreting neutron stars 
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ABSTRACT 

We examine the hydro-magnetic stability of magnetically confined mountains, which 
arise when material accumulates at the magnetic poles of an accreting neutron star. 
We extend a previous axisymmetric stability analysis by performing three-dimensional 
simulations using the ideal-magnctohydrodynamic (ideal-MHD) code ZEUS-MP, inves- 
tigating the role played by boundary conditions, accreted mass, stellar curvature, and 
(briefly) toroidal magnetic field strength. We find that axisymmetric equilibria are 
susceptible to the undular sub-mode of the Parker instability but are not disrupted. 
The line-tying boundary condition at the stellar surface is crucial in stabilizing the 
mountain. The nonlinear three-dimensional saturation state of the instability is char- 
acterized by a small degree of nonaxisymmetry (< 0.1 per cent) and a mass ellipticity 
of e ~ 10~ 5 for an accreted mass of M a = 10~ M . Hence there is a good prospect 
of detecting gravitational waves from accreting millisecond pulsars with long-baseline 
interferometers such as Advanced LIGO. We also investigate the ideal-MHD spectrum 
of the system, finding that long-wavelength poloidal modes are suppressed in favour 
of toroidal modes in the nonaxisymmctric saturation state. 
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1 INTRODUCTION 

There exists strong observational evidence that the mag- 
netic dipole moment of accreting neutron stars in X- 
ray binaries, /i, decreases with accreted mass, M a 



(|Taam & van de Heuvell 


1986: van den Heuvel & Bitzarakil 


1995!) , although Wiiers 


|l997|) noted that extra vari- 



ables may enter this relation 
for field reduction have been 
cclcratcd 



Numerous mechanisms 
proposed, such 



Ohmic decay 



a s ac- 

Konar fc Bhattacharval 1 19971 : 
lUrpin fc Konenkovl ll997T ). vortex-fluxoid interactions in 
the superconducting core (|Muslimov fc Tsyganl 1 19851 : 



ISrinivasan et al. 1990t). and magnet i c screening or burial 
llBisnovatvi-Kogan fc Komberel 1974 : Romanilll990l: IZhand 



Il998l : IPavne fc Melatosl 1200 



Lovelace et aU I2005T ). The 



reader is referred to iMelatos fc Phinnev (120011 ) for a com- 
parative review. 

Magnetic burial occurs when accreted plasma flowing 
inside the Alfven radius is chanelled onto the magnetic poles 
of the neutron star. The hydrostatic pressure at the base 
of the accreted column overcomes the magnetic tension of 
the magnetic field lines and spreads eq uatorwards, thereby 
distorting the frozen-in magnetic flux (|Melatos fc Phinnevl 



E-mail: mvigeliu@physics.unimelb.edu.au 



l200ll ). |Pavne fc Melatosl (|2004h . hereafter PM04, computed 
self-consistently the unique quasistatic sequence of ideal- 
MHD equilibria that describes how burial proceeds as a 
function of M a , while respecting the flux freezing constraint 
of ideal MHD. They found that the magnetic field is com- 
pressed into an equatorial belt, which confines the accreted 
mountain at the poles. A key result is that \i is reduced sig- 
nificantly once M a exceeds ~ 10~ 5 Mp), five orders of mag- 
nitude above prev ious estimates jBrown fc Bildstenl 1 19981 : 
iLitwin et al"1l200ll ). 

Generally speaking, one expects highly distorted hydro- 
magnetic equilibria like those in PM04 to be disrupted on 
the Alfven time-scale by a plethora of MHD i nstabilities. 
Surprisingly, however, IPavne fc Melatosl (|2007l ) (hereafter 
PM07) found the equilibria to be stable to axisymmetric 
ideal-MHD modes. When kicked, the mountain performs ra- 
dial and lateral oscillations corresponding to global Alfven 
and compressional modes, but it remains intact. 

An axisymmetric analysis, however, excludes insta- 
bilities involving toroidal modes and is therefore incom- 
plete. In general, three-dimensional effects alter MHD 
stability, quantitatively and qualitatively. For example, 
iMatsumoto fc Shibatal (| 19921 ) found the growth rate of the 
three-dimensional Parker instability t o be higher th a n that 
of its two-dimensional counterpart. iMasada et al.l (|2006l ) 
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proved that newly born neutron stars containing a toroidal 
field are stable to the axisymmetric magneto-rotational in- 
stability yet unstable to its nonaxisymmetric counterpart. 
Differences between the two- and three-dimensional stabil- 
ity of MHD equilibria ar e also observed in a variety of so- 
lar contexts (|Priestlll984h and in tokamaks (|Lifschit zl 1 19891 ; 
iGoedbloed fc Poedtsll20o3 ). 

The central aim of this paper is to perform fully three- 
dimensional, ideal-MHD simulations to assess the stability 
of magnetic mountains, generalizing PM07. Importantly, we 
compute not just the linear growth rate but also the non- 
linear saturation state of any unstable modes. The latter 
property is what matters over the long accretion time-scale 
when evaluating magnetic bur ial as the cause of the observed 
reduction in p (|Pavnel [2005f ). the persis tence of millisec- 
ond o scillations in type-I X-ray bursts (|Pavne fc Melatosl 
l2006bl ). and gravitational radiation from magnetic m oun- 



tains (jMelatos fc Pavnell2005l ; iPavne fc Melatosl 1200681 ) . 

The structure of the paper is as follows. We introduce 
our numerical setup in section [2] and validate it against pre- 
vious axisymmetric results in section [3] characterizing the 
controlling influence of the boundary conditions for the first 
time. In section[4] we present three-dimensional simulations, 
which display growth of unstable toroidal modes. The insta- 
bility is classified according to its dispersive properties and 
energetics, and the nonlinear saturation state is computed as 
a function of M a . We compute the spectrum of global MHD 
oscillations in section [S] Resistive effects are postponed to a 
future paper. 



2 NUMERICAL MODEL 

The accretion problem contains two fundamentally differ- 
ent time-scales: the long accretion time (~ 10 s yr) and the 
short Alfven time (~ 10 -3 s). The wide discrepancy pre- 
vents us from treating the accretion problem dynamically, 
i.e. in a full MHD simulation, where mass is added through 
the outer boundary onto an initially dipolar field. Instead, 
for a given value of M a , we compute the magnetohydro- 
static equilibria, using the Grad-Shafranov solver developed 
by PM04, then lo ad it into the ideal-MHD solver ZEUS-MP 
jHaves et al.ll2006l ) to test its hydromagnetic stability on the 
Alfven time-scale ta- We find below that all quantities reach 
their saturation values after ~ Wta ~ 10~ 2 s at a particu- 
lar value of M a (e.g. 10 _4 Mq). As M a changes slowly, over 
~ 10 s yr, the saturation values adjust in a quasistatic way 
on the Alfven time-scale. In practice, to study a different 
value of M a numerically, we recalculate the Grad-Shafranov 
eqilibrium and load the new equilibrium into zeus-mp. 



2.1 Magnetic mountain equilibria 

Analytic and numerical recipes for calculating self-consistent 
ideal-MHD equilibria for magnetic mountains are set out 
in PM04. Here, we briefly restate the main points for the 
convenience of the reader. 

An axisymmetric equilibrium is generated by a scalar 
flux function ip(r,6), such that the magnetic field 



automatically satisfies V-B = 0. We employ the usual spher- 
ical coordinates (r,9,4>), where 8 — corresponds to the 
symmetry axis of the magnetic field before accretion. In 
the static limit, the mass conservation and MHD induction 
equations are identically satisfied, while the component of 
the momentum equation transverse to B reduces to a sec- 
ond order, nonlinear, elliptic partial differential equation for 
ip, the Grad-Shafranov (GS) equation (PM04): 



A 2 ip 
with 



-F'(ip) exp[-(v? - v?o)/c 2 ], 



por 2 sin 2 I 



dr' 



+ 



t) 



<:) 



08 \sin8d8 



(2) 



(3) 



Formally, F(ip) is an arbitrary function. However, the 
ideal-MHD flux- freezing constraint, that matter cannot 
cross flux surfaces, imposes an additional conservation law 
on the mass-flux ratio AM /dtp, 



cj AM 
2^"dvT 



Asr sm8\S7ip\ x e 



-(<p-<f>o)/ c a 



(4) 



which determines F(ip) uniquely when solved simultane- 
ously with @. The integration in @ is performed along 
the field line ip — const, <p = const. The solution is insen- 
sitive to the exact form of AM /dip; we thus distribute the 
accreted mass M a uniformly over ^ ip ^ ip a , where ip a 
is the flux enclosed within the polar cap, and take ip to be 
dipolar, initially, with hemispheric flux ip*. 

In writing (0 and © , we approximate the gravitational 
field as uniform over the height of the mountain, and hence 
write the gravitational potential tp as 



<fi = GM«r/Rt 



(5) 



with ipo = GM*/R*, where M* and R* denote the stellar 
mass and radius respectively. We also assume an isothermal 
equation of state, p = c 2 p, where c s denotes the sound speed. 

By wor king in the ideal-MHD limit, we neglect elas- 
tic stresses (jMelatos fc Phinnevl l200ll; lHaskell et al.l 120061 ; 
OwenlbOOrih the Hal l drift (iGeppert fc Rheinhardtl [2051 ; 
Cumming et all 120041; IPons fc Geppertl 120071), and Ohmic 
diffusion (|Romanilll99d ; IGeppert fc Urpinlll99i ). In partic- 
ular, Ohmic diffusion causes the mass quadrupole moment 
of the magnetic mountain (and p) to saturate above a cer- 
tain value of M a and may also affect the stability to resistive 
MHD (e.g. ballooning) modesQ We defer investigating these 
resistive effects to a forthcoming paper [Vigelius & Melatos 
(in preparation)]. 

We solve and @ simultaneously using the relax- 
ation algorithm described in PM04, subject to the bound- 
ary conditions ip(R*,9) — ip* sin 2 8 (line tying at the sur- 
face), dip/dr(R m ,0) = at the outer boundary r = R m , 
ip(r,0) = 0, and dip/d8(r,n/2) = (north-south symme- 
try). 



2.2 Evolution in ZEUS-MP 

In this paper, we explore numerically how the axisymmetric 
GS equilibria evolve when subjected to a variety of initial 



(1) 



1 J. Arons, private communication. 
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and boundary conditions in three dimensions. To achieve 
this, we employ the parallelized, gene ral purpose, time - 
dependent, ideal- MHD solver ZEUS-MP (|Haves et al]|2006h . 
ZEUS-mp integrates the equations of ideal MHD, discretized 
on a fixed staggered grid. The hydrodynamic part is based 
on a finite-difference advection scheme accurate to second 
order in time and space. The magnetic tension force and the 
induction equation are solved via the meth od of characteris- 
tics a nd constrained transport (MOCCT) (|HawIev fe Stonel 
Il995t ). whose nume rical implementatio n in zeus-mp is de- 
scribed in detail bv lHaves et al.l (|2006T ). 

We initialize ZEUS-mp with an equilibrium computed 
by the GS code and described by B(r, 8) and p(r,9), ro- 
tated about the z axis to generate cylindrical symmetry. We 
introduce initial perturbations by taking advantage of the 
numerical noise produced by the transition between grids in 
the GS code and zeus-mp. 

We adopt dimensionless variables in zeus-mp satisfying 
fio = G = c s = ho = 1, where ho = c 2 R 2 /GM, denotes the 
hydrostatic scale height. The basic units of mass, magnetic 
field, and time are then M = h c 2 JG, B = [c 4 /(G7io)] 1/2 , 
and ro = ho/c s . The grid and boundary conditions are spec- 
ified in appendix 1X1 



2.3 Curvature rescaling 

In general, the characteristic length-scale for radial gradi- 
ents (ho) is much smaller than the length-scale for latitu- 
dinal gradients i?„, creating numerical difficulties. However, 
in the small- M a limit, it can be shown analytically (PM04, 
PM07) that the structure of the magnetic mountain depends 
on Rt, and M* through the combination ho oc i? 2 /M*, not 
separately. We therefore artificially reduce R, and M* , while 
keeping ho fixed, to render the problem tractable computa- 
tionally. It is important to bear in mind that invariance of 
the equilibrium structure under this curvature rescaling does 
not imply invariance of the dynamical behaviour, nor is it 
necessarily applicable at large M a . 

A standard neutron star has M* = 1.4M©, R* = 10 6 
cm, B* = I0 12 G, and c s = I0 8 cm s" 1 , giving ho = 53.82 
cm, a = R*/ho = 1 .9x I0 4 , and ro = 5.4x ICT 7 s. We rescale 
the star to M» = f.O x I0" 5 M Q and R', = 2.7 x I0 3 cm, 
reducing a to 50 while keeping it large. The base units for 
this rescaled star (see section |2~2")) are then Mo — 8.1 x I0 24 
g, po = 5.2 x I0 19 g cm" 3 , Bo = 7.2 x I0 17 G, and r = 

5.4 x 10~ 7 s. The critical accreted mass above which the 
star's magnetic moment starts to change, M c , is defined by 
equation (30) of PM04: 

A characteristic time-scale for the MHD response of the 
mountain is the Alfven pole-equator crossing time, ta = 
7tJ?*/(2wa), where va = (B 2 / pop) 1 ^ 2 is the Alfven speed. 
Clearly, va is a function of position and time, so the defini- 
tion of ta is somewhat arbitrary. Typically, at the equator, 
we find B ~ IO _6 -Bo and p ~ 10 -11 po, empirically implying 
r A « 250r . 



Table 1. Simulation parameters. M a is the accreted mass, in 
units of the characteristic mass M c (section 12.31 1 , and a = R„/ho 
measures the curvature of the rescaled star (section 12. 3D . The 
conditions at the outer boundary (r = R m ) are either outflow 
(zero gradient in all field variables) or inflow (pinned magnetic 
field); cf. also appendix lAl 



Model 


M a /M c 


a 


Axisymmetry 


Boundary 


A 


1.0 


50 


yes 


outflow 


B 


1.0 


50 


yes 


inflow 


D 


1.0 


50 


no 


outflow 


B 


1.0 


50 


no 


inflow 


F 


0.6 


50 


no 


outflow 


G 


1.4 


50 


no 


outflow 


J 


1.0 


75 


no 


outflow 


K 


1.0 


100 


no 


outflow 



3 AXISYMMETRIC STABILITY AND 
GLOBAL OSCILLATIONS 

PM07 demonstrated the axisymmetric stability of magnetic 
mountains using the serial ideal-MHD solver ZEUS-3D. Here, 
we start by repeating these axisymmetric simulations in the 
parallel solver zeus-mp, in order to verify the mountain im- 
plementation in ZEUS-3D and zeus-mp, generate an axisym- 
metric reference model, and understand the effect of the 
boundary conditions, which were not investigated fully in 
previous work. The simulation parameters are detailed in 
Tabled] (models A and B). 

3.1 Reference model 

Model A, in which we set M a /M c = f.O and the outer 
boundary condition to outflow, serves as a reference case. 
Fig.[T]displays a time series of six r-9 sections for ^ t/rA ^ 
3.6, showing density contours (dashed curves) and the mag- 
netic field lines projected into the plane = (solid curves). 
The axisymmetric equilibrium (top-left panel) reveals how 
the bulk matter is contained at the magnetic pole by the 
tension of the distorted magnetic field. 

The mountain in Fig. [1] performs damped lateral oscil- 
lations without being disrupted. The (unexpected) stability 
of this configuration is due to two factors. First, the config- 
uration is already the final, saturated state of the nonlin- 
ear Parker instability, whi ch is reached quas istatically dur- 
ing slow accretion (PM04; [Mouschoviasll 19741 ). Second, line- 
tying of the magnetic field at r = R* significantly changes 
the structure of t he MHD wave spectrum i n a wa y that 
enhances stability. iGoedbloed fc Halberstadtl (|l994T ) found 
that, in a homogenous plasma, a superposition of Alfven- 
and magnetosonic waves is needed to satisfy the line-tying 
boundary conditions. As a consequence, the basic, unmixed 
MHD modes are not eigenfunctions of the linear force oper- 
ator, and thus the spectrum is modified. 

The mass quadrupole moments, 

Qij = J d 3 x (Zx'iXj - r' 2 %)p(x'), (7) 

of the mountain in Fig. [1] are plotted versus time in Fig. [2] 
We note first that Q\2 = and Q22 = — Q33/2, as expected 
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52 54 56 58 52 54 56 58 52 54 56 58 

r/h r/h r/h 



Figure 1. Meridional section of model A at 4/ta = 0,0.4,1,1.8,2.4,3.6 (top left to bottom right). Shown are density con- 
tours (dashed curves) with values log 10 (p/pQ) = —13, —12, —11, —10.7, —10.5, —10.3, and flux surfaces with footpoints at r = Ft*, 
6 = 0.10, 0.12, 0.15, 0.20, 0.39, 0.79 (solid curves). Lateral oscillations of the equatorial field lines are clearly visible. The outflow bound- 
ary condition at r = R m makes the field lines flare towards the magnetic pole. This is visible most clearly for the line whose footpoint 
lies at = 0.39. 



for an axisymmetric system. We can comp are Fig. [2]directly 
with th e ellipticity e oc Q33 computed bv lPavne fc Melatosl 
(2006a) Q. These authors found two dominant global modes, 
an Alfven and an acoustic mode, claimed to be analogous to 
the fundamental modes of a gravitating, magnetized plasma 
slab. We cannot resolve the compressional modes in Fig. [2] 
but the latitudinal (Alfven) mode is clearly visible through 
the oscillations in Q22 and Q33. 

An oscillation cycle proceeds as follows. The first mini- 
mum of Q33, and hence e, at t — 0.4ta in Fig.[2j corresponds 
to the top right panel of Fig. [T] The mountain withdraws ra- 
dially and poleward. Polar field lines move closer to the mag- 
netic pole, while equatorial field lines are drawn towards the 
equator. At r = ta, the mountain spreads and Q33 reaches 
a maximum in Fig. [2] The damping observed in Fig. [2] arises 
solely from numerical dissipation; neither viscosity nor re- 
sistivity are included in our version of zeus-mp. 

3.2 Outer boundary 

The "flaring up" of magnetic field lines at the pole, observed 
by PM07, is an artifact of the outflow boundary condition 

2 IPavne fc Melatosl l|2006ah simulated a polar cap with b = 10 as 
against b = 3 in model A. 




Figure 2. Mass quadrupolc moments for model A, normalised to 
the maximum of Q33 (1.3 X 10 25 g cm 2 ), as a function of time 
in units Alfven time. We find Q22 = — Q33/2 and Q12 = 0, as 
expected for an axisymmetric configuration. The global Alfven 
oscillation is damped by numerical viscosity. 

at r = R m . In order to check this, we repeat the simulation 
of model A but switch to an inflow boundary condition 
(model B). The density distribution is similar in the two 
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52 54 56 58 52 54 56 58 

r/h r/h B 



Figure 3. Meridional section of model B, showing density (left) 
and projected magnetic field (right) at t/r\ = 3.6 (solid curves). 
The field lines are pinned to the outer r boundary by the inflow 
boundary condition. For comparison, the density and magnetic 
field of model A are overplotted (dotted curves). 

0.4 f ' 1 1 1 : 

0.2 r i 



& 0.0 - 
-0.2 7 




Figure 4. Mass quadrupolc moments for model B, normalised 
to the maximum of Q33 (1.3 X 10 25 g cm 2 ) as a function of time 
in units of the Alfven time. We find Q22 = — Q33/2 and Q12 = 
0, as expected for an axisymmetric configuration. The mountain 
performs damped lateral oscillations with twice the frequency of 
model A. 

models, as is clear from Fig. [3] However, the inflow BC 
artifically pins the magnetic field to the outer boundary, 
introducing magnetic field discrepancies (mostly in the outer 
layers, where p is negligible). 

Alfven waves are similar to transverse waves on a string, 
where the magnetic tension provides the restoring force. 
Models A and B are therefore equivalent to a vibrating string 
with one end free and fixed respectively. We expect the os- 
cillation frequency of the fundamental mode in model B to 
be twice that of model A. This is indeed observed in the 
oscillations of the quadrupole moments, displayed in Fig. [4] 
Q22 and Q33 in Fig. |4]oscillate at 0.6 times the period in Fig. 
[2] A comprehensive analytic computation of the MHD spec- 
trum, including discrete and continuous components, will be 
attempted in a forthcoming paper. 

3.3 Uniform toroidal field 

There are strong theo retical indications that the magne- 
torotational instability (|Balbus fc Hawle"vlll99cf ) acts during 




Figure 5. Model A repeated with the same parameters as in Ta- 
bic [l] including a uniform toroidal field = 10 -7 So = 0.35B P . 
The mountain is defined by the orange isosurface p(r, 8, <fr) = 
1.04 X 10 9 g cm 3 , while red denotes the neutron star surface 
r = Ft*. In order to improve visibility, all length scales of the 
mountain and the field lines (blue) are magnified five-fold. The 
field exhibits a helical topology, which is most distinct in the polar 
flux tubes where the poloidal contribution is weakest. 



core collapse supernova explosions to generate a substan- 
tial toro idal field component B& ~ B v be neath the stellar 
surface (|Cutlerll20o3 ; lAkivama et al.ll200ot ) . The hydromag- 
netic stability of equilibria with 7^ will be discussed 
thoroughly in a forthcoming paper. In this subsection, for 
completeness, we present the results of a preliminary inves- 
tigation. 

Let us rerun model A with the same initial conditions 
while applying a uniform B$ — 1Q~ 7 Bq = 0.35B P through- 
out the integration volume and at r — R, , where the poloidal 
field component is defined as B p = (B 2 +B|) 1//2 , taken at the 
point [x= (r-R*)/h ,9,(j>] = (1(T 3 , 0.012, 0). B^ is uniform 
only initially and is allowed to evolve nonuniformly as ZEUS- 
mp proceeds. This procedure leads to a non-equilibrium con- 
figuration, because we do not generalise and solve again the 
GS equation @ to accomodate B^ 7^ 0. Nevertheless, it pro- 
vides us with some insight into the stability of a field with 
nonzero pitch angle. 

Fig. [S] displays the result of this numerical experiment 
after t = 7 At a- The toroidal field component creates a he- 
lical field topology in the polar region far from the surface, 
where the poloidal field is comparably weak. In the equa- 
torial region, however, the poloidal field is still dominant 
and the structure remains unchanged from Fig. [1] Remark- 
ably, the toroidal field does not alter the stability of the 
system qualitatively, at least for the parameters of model A 
[We expect a st r onger effect in other p a ramet er regimes; see 
iLifschitj (|l989T l; lGoedbloed fc Poedtj (|2004T l]. The elliptic- 
ity, displayed in Fig. [6] (solid curve), exhibits characteristic 
oscillations with a period ~ 30 per cent smaller than that of 
the purely poloidal configuration (dotted curve), which can 
be explained simply by the increase in the Alfven speed. In 
addition, the saturation ellipticity is ~ 3 per cent higher 
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t/r, 




Figure 6. Mass ellipticity of model A with (solid) and without 
(dotted) a uniform B^. The toroidal field component leads to a 
shorter oscillation period and to a higher saturation ellipticity. 



Figure 8. Evolution of the azimuthal magnetic field component 
B^ at (x , B) = (1(T 3 , 0.01) in model D as a function of longitude 
<j> (in radians) and time t (in units of the Alfven time.) 



than in model A; the magnetic tension increases with B, 
sustaining the mountain at a lower colatitude. 



4 NONAXISYMMETRIC STABILITY 

We turn now to the three-dimensional evolution of a mag- 
netised mountain in zeus-mp. The chief finding, presented 
below, is that the initial (axisymmetric) configuration be- 
comes unstable to toroidal perturbations, but that, after a 
brief transition phase, the system settles into a new (nearly 
axisymmetric) state, which is stable in the long term. Sec- 
tion [4TT] compares the results to the axisymmetric reference 
model A. The magnetic and mass multipole moments are 
computed in section |4~21 the influence of the boundary con- 
ditions is considered in section 1131 and the component-wise 
evolution of the energy is examined in section 14.41 A scaling 
of the mass quadrupole moment versus M a is derived em- 
pirically in section l4~5l The curvature rescaling is verified in 
section 14.61 



4.1 General features 

Model D starts from the same configuration as model A 
(M a /M c — 1.0, outflow at r = R m ) but is evolved in three 
dimensions. Six snapshots of a density isosurface (orange) 
and magnetic field lines (blue) are depicted in Fig. [7] At 
r ~ 0.8ta = 200to, the system undergoes a violent tran- 
sition. The field lines bend in the <j) direction, indicating 
that the initial axisymmetric configuration is unstable to 
toroidal modes, a channel that is evidently not present in 
axisymmetric simulations. This hypothesis is supported by 
Fig. [5] which plots B$ at x = 10~ 3 and = 0.01 as a 
function of cj> and t. The magnetic field takes the form of 
an azimuthal travelling wave B^, ~ exp[i(m0 — ut)]. From 
Fig. [S] we measure the phase speed to be approximately 
v p = ujR- t /m — 27 'v a, where the Alfven speed va is mea- 
sured at (x,0,4>) = (10 -3 , 0.01, 0.1) and we assume m = 1. 
An inhomogenous plasma generally supports mixed magne- 
tosonic/ Alfven modes, so ity does not necessarily equal v\ 
or the fast /slow magnetosonic speed. The magnitude of B$ 




Figure 9. \B^\ at (x,d,cj>) = (10 -3 , 1.5, 1.3) for model D, simu- 
lated at higher resolution (JV^ = 32, solid) and lower resolution 
(JV^ = 8, dashed). Small wavelength perturbations grow faster. 
We track the absolute value of | B^ \ in order to isolate better the 
dominant mode, as the system exists in a superposition of stable 
and unstable modes, and Ba, switches sign. 



is comparable to the magnitude of the polar magnetic field 
B p = 2.9 x 10~ 7 B . 

The deviations of the mountain isosurface, defined by 
p(r,8,(f)) = 1.04 x 10 9 g cm 3 , from axisymmetry are small 
(< 5 per cent laterally and < 0.06 per cent radially during 
the transition phase). The isosurface spreads outward by 
~ 32 per cent relative to ites initial position. 

Fig. [5] demonstrates how the instability grows. We plot 
\B^\ at the (arbitrary) position (r,6,(j>) = (R,, 1.5, 1.3) ver- 
sus time. The solid curve corresponds to a higher toroidal 
resolution (Nj, = 32 grid cells in (f> direction) than the 
dashed curve (N$ = 8). We note first that \B^\ grows expo- 
nentially with time, as expected in the linear regime. The 
growth rate is measured to be F = Im(ai) = 0.05^ 1 — 
12.5r^ 1 . Second, the instability is manifestly associated with 
toroidal modes. The magnetic perturbation, <5B, induced by 
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Figure 7. Density and magnetic field of model D at t/t\ = 0, 1,2,3,4,5 (from top left to bottom right). The mountain is defined by 
the orange isosurface p(r, 9, <j>) = 1.04 X 10 9 g cm 3 , while red indicates the neutron star surface r = R*. In order to improve visibility, all 
length scales of the mountain and the field lines are magnified five- fold. The mountain becomes unstable to toroidal modes at r fss 0.8ta- 
It subsequently relaxes to a new nonaxisymmetric equilibrium. The footpoint of the blue fieldlines is at the stellar surface while green 
ficldlines are traced starting from the equator. Green field lines eventually become topologically disconnected (see text). 



a linear Lagrangian displacement £ is <5B = V x (£ x B). By 
writing out the vector components, one sees that 8B4, =^ 
implies £^ 7^ 0, provided the unperturbed field has the form 
B = B r {r,e)e r + Be(r,8)e g . 

The dashed curve in Fig. [9] tracks \B$\ for a simulation 
carried out at a lower resolution (N<p = 8). The instabil- 
ity grows significantly slower with F — Im(oj) = 0.02r _1 = 
5r^ . We conclude that F scales with the wavelength A of 
the perturbation roughly as A -1 / 2 . The piecewise-straight 
appearance of the dashed curve in Fig. [9] shows that the 
global oscillations are governed by a superposition of unsta- 
ble (growing) and stable wave modes. 

What type of instability is at work here? In order to an- 
swer that question, we first write down the change in poten- 
tial energy associated with a Lagrangian displacement £ in a 
form that reveals the physical meaning of the different con- 
tribu tions (Biskamp 1993; lLifschit jfl989l ; lGreene fc Johnson! 
1 19681 ): 



8W V 



dV [\Q±\ 2 + \V ■ £± + 2k ■ £ A 



\ 2 Bl 



+clp\V ■ i\ 2 - j^Cx xb)-Q 

-2(f ± • Vp)(« • Cx) - (£* • V^)V • (pO] , (8) 



We i nclude the term due to gravity (|Goedbloed fc Poedtsl 
120041 ) and define j = (ij'V x B (current density), Q = 
V x (£ x B) (change in B as a response to £), n = (b ■ V)b 
(field line curvature), and b = B/_B. The subscripts _L, || 
refer to the magnetic field, such that a± = a — (a ■ b)b 
and a || = a • b. The first three (stabilising) terms are the 
potential energy of the shear Alfven mode, the fast mag- 
netosonic mode, and the (unmagnetized) sound mode. They 
are all positive definite. The last three terms may have either 



sign. The term proportional to j« causes the current-driven 
instabilities, the curvature term causes pressure-driven in- 
stabilities (when k ■ V p > 0), and the final term causes 
gravitational instabilities. 

We first note that ju = 0, ruling out current-driven in- 
stabilities. Furthermore, we have k = K r (r, 0)e r + ne(r, 9)ee 
in our particular field geometry. In principle, this term ad- 
mits pressure-driven instabilities. However, we would expect 
such instabilities, if they exist, to also grow in an axisymmet- 
ric system, yet they do not. This suggests that the instability 
we see in Figs. [7}{9]is associated with a toroidal dependence 
in £, leaving the gravitational term, which indeed contains 
c^,£ contributions. 

One prominent gravitational mo de is the Parker o r 
magnetic buoyancy instability (PM04; 



Its physics was elucidated by iHughes 



Mouschovias 197< 
Cattaneol (Il98 



D 

for a plane-parallel, stratified atmosphere with a horizon- 
tal field increasing with depth z. The instability involves 
an interchange sub-mode and an undular sub-mode. The 
interchange sub-mode satisfies k y = £ y = 0. We do not ob- 
serve this mode in our system because (i) it should also 
be present in two dimensions, as it does not rely on a 
toroidal dependence, yet it is absent; and (ii) it is incon- 
sistent with the line-tying boundary condition at r = _R*. 
On the other hand, undular modes compress the plasma 
along field lines, even in systems which are interchange sta- 
ble. In two dimensions, they are restricted to k x = Q x = 0, 
whereas a non-vanis h ing Q x is allowed in three dimensions. 
IHughes fc Cattaneol (|l987f ) showed that SW P is minimized 
for k x — > 00, consistent with the results in Fig. [9] the insta- 
bility grows faster, if we allow smaller wavelength perturba- 
tions by increasing iV^. 

When va is uniform the growth rate of the Parker in- 
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stability reaches an asymptotic maximum Fp ~ (g/A) 1 / 2 for 
k x A,k y A 2> 1. Here, A = v%/g is the scale height for a strat- 
ified atmosphere with uniform gravitational acceleration g. 
We recognize (A/g) 1 ' 2 as the characteristic free fall time over 
one scale height. In the units specified in section ROl we find 
TpTo — A — g — 1, two orders of magnitude higher than the 
observed growth rate F ~ lO -2 ^ 1 . The discrepancy arises 
because the Parker instability cannot grow freely in the belt 
region, since the adjacent plasma at higher latitudes effec- 
tively acts as a line-tying boundary for the magnetic field. 

The snapshot at t = ta (top-middle panel in Fig. [7J) 
demonstrates how the instability starts in the equatorial re- 
gion, whose magnetic belt represents the endpoint of the 
two-dimensional Parker instability. The undular Parker sub- 
mode releases gravitational energy by radial plasma flow 
towards the neutron star's surface. At the same time, the 
magnetic field is rearranged such as to minimize the (radial) 
gradient in B. Importantly, the undular mode is not avail- 
able in the axisymmetric case. The extra degree of freedom 
in the (j> direction allows perturbations to develop which do 
no work against the magnetic pressure, destabilising the belt 
region. 

Particularly interesting here is the formation of topolog- 
ically disconnected field lines (green curves in Fig.0. These 
occur when field lines are pushed out of the radial boundary 
surface. They are then disrupted and can subsequently re- 
connect at the equatorial boundary, forming O-type neutral 
points ( "bubbles") and associated Y-type points. It is impor- 
tant to note, however, that the formation of these bubbles 
is not an unphysical boundary effect. Instead, in a realistic 
setting, even a small resistivity leads to reconnection and 
thus to a topological rearrangement o f the field. The effect i s 
similar to the formation of plasmoids (|Schindler et al.ll 19881 ) . 
B r switches sign at the magnetic equator, implying the ex- 
istence of a current sheet. Reconnection then leads to the 
creation of magnetic X-type neutral points and the associ- 
ated bubbles. We discuss resistive effects and the importance 
of these bubbles to resistive instabilities in an accompanying 
paper. 

The concept of rational magnetic surfaces, where the 
field lines close upon themselves, plays an important role in 
a local plasma stability analysis. The bending of field lines 
as a result of a Lagrangian displacement is associated with 
an increase in potential energy. Hence, for almost all insta- 
bilities to occur, this contribution, which can be expressed 
as B ■ V£, needs to be small. In a tokamak geometry, it can 
be shown that this term vanishes on a rational surface. The 
spatial location of rational surfaces is directly related to the 
pitch angle B^j B p . We defer a detailed analysis of the ratio- 
nal magnetic surfaces in our problem to a forthcoming paper 
and restrict ourselves to a brief discussion in the following 
paragraph. 

In Fig. 1101 we plot the pitch angle as a function of the 
coordinate r\ along the field line (right panel) for four differ- 
ent field lines (left panel) in model D. We first note that the 
toroidal component stays below 20 per cent of the poloidal 
component along all four field lines. Furthermore, the ab- 
solute magnitude of the pitch angle tends to increase with 
colatitude. This is consistent with the previous discussion. 
The Parker instability (and hence B$) dominates close to 
the magnetic equator. The wave-like character of the insta- 
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Figure 11. Mass quadrupolc moments for model D, normalised 
to the maximum of Q33 (1.30x 10 25 g cm 2 ) as a function of time in 
units of the Alfven time. The sytem develops a substantial asym- 
metry, characterised by the off-diagonal element Q12, during the 
relaxation phase, before settling down to a nearly axisymmetric 
state. 



bility is vividly demonstrated by the zero crossings of the 
pitch angle. 



4.2 Mass and magnetic quadrupole moments 

At t » 2ta, the system in Fig. [7] settles down to a stable 
state which differs from the initial configuration, primarily 
by being nonaxisymmetric with respect to the pre-accretion 
magnetic axis. The field lines whose footpoints are at a low 
colatitude move towards the magnetic poles. This behaviour 
is reflected in the mass quadrupole moments, plotted against 
time in Fig. [11] The transition to a nonaxisymmetric mag- 
netic field configuration at t ~ ta is accompanied by a sud- 
den rise in the off-diagonal moment Q12. However, by the 
time the mountain settles down at t ~ 2ta, axisymmetry is 
largely restored and Q12 decreases. Fig. [TT] shows that Q12 
oscillates before damping down, with a remarkably low de- 
viation from axisymmetry of Q12/Q33 < 0.1 per cent in the 
final state. 

We reiterate that the final state is not the same as the 
initial state, even though it is nearly axisymmetric. Further- 
more, the final state is stable. This is the main result of the 
paper, as far as astrophysical applications are concerned. 

Why does Qij decrease? Naively, one would not expect 
a significant change, given that the Parker instability pre- 
dominantly acts in the equatorial belt region, while most 
of the plasma is located at the magnetic pole. The answer 
can be found in the Lorentz force, which balances the lateral 
pressure gradient. Fig. [T2]shows how the relative strength of 
the toroidal component of the Lorentz force, (J x B)^/| J x B| 
(dashed curve), grows as a function of time in model D. As 
B^ grows, following the onset of the instability, the force per 
unit volume develops a toroidal component while its lateral 
component decreases. Hence the (approximately unchanged) 
lateral hydrostatic pressure gradient forces the mountain to 
slip towards the equator. After the system settles down, B^ 
decreases and the lateral components of J x B and Vp read- 
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Figure 10. Magnetic pitch angle B^/Bp (right panel) as a function of the coordinate rj along four magnetic field lines 1-4 for model D, 
for a snapshot taken at t = 5ta ■ The positions of the field lines are depicted in the left panel. The mountain is defined by the orange 
isosurface p(r, 0,<f>) = 1.04 X 10 9 g cm 3 , while red denotes the neutron star surface r = R» . In order to improve visibility, all length scales 
of the mountain and the field lines (blue) are magnified five-fold. 




Figure 12. Meridional section for model D at t/rA = 0, 1, 2, 3, 4, 5 (top left to bottom right). Shown are the density contours (dashed 
curves) for log 10 (p/p' ) = —13, —12, —11, —10.7, —10.5, —10.3, and the normalised Lorentz force per unit volume (J X B)^/|J X B| 
(solid curves) for the values 0.1, 0.5, 0.9. The Lorentz force develops a toroidal component as B^ increases, but its poloidal component 
diminishes, allowing the poloidal pressure gradient to push the mountain equatorwards. 



just to balance each other, leading to the stable equilibrium 
state. 

The ellipticity e oc Q22 reaches a local maximum during 
the transition phase at t ~ ta and subsequently drops. The 
mass quadrupole moment of the final configuration is ~ 33 
per cent lower than in model A. The asymptotic values of 
Qij for the eight models in Table [TJ are tabulated in Table 
[2] normalized to Qsz(t = 8ta) for model A. 



The nonvanishing components of the magnetic dipole 
and quadrupole moments di m (r = R m ), defined as 

di m {r = Rm) = R\t 2 J dfi YC m B r (9) 

(see appendix [Cj , are displayed as functions of time in Fig. 
1111 The dipole moment dio = 4(7r/3) 1//2 /i increases rapidly 
during the transition phase, reaching an asymptotic maxi- 
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Table 2. Asymptotic values of Qij for the eight models in Table 
[T] normalized to Q33 = Q33(t = 8ta) for model A. We select 
t = 8r A (models A-E), t = 5t a (models F & G), and t = 4t a 
(models J & K) to compute the asymptotic value. 



Model 


Q12/Q33 


Q22/Q33 


Q33/Q33 


A 


-3.1 x 10" 8 


-0.50 


1.00 


13 


-3.1 x 10" 8 


-0.51 


1.02 


D 


1.6 x 10" 3 


-0.21 


0.64 


E 


-3.1 x icr 8 


-0.51 


1.02 


F 


-1.9 x 10" 3 


-0.23 


0.60 


G 


8.7 x icr 4 


-0.17 


0.68 


J 


5.7 x 10" 2 


-1.27 


3.9 


K 


2.5 x 10" 4 


-4.72 


12 



Table 3. Non- vanishing components of the asymptotic mag- 
netic dipole moments dio /R 3 ^ and magnetic quadrupole moments 
qw/Rm, both normalised to dio/R^ = dio(t = 8rA)/i? 3 n . We se- 
lect t = 8ta (models A-E), t = 5ta (models F & G), and t = 4ta 
(models J & K) to compute the asymptotic value. 



Model 


dw/d\o 


Re(d 21 )/(d 10 R m ) 


lm(d 2 i)/(dioR m ) 


A 


1.0 


-4.6 x 10" 8 


-4.8 X 10- 9 


U 


1.3 


-5.9 x 10~ 8 


-6.1 x 10~ 9 


D 


7.4 


4.8 x 10" 2 


-4.3 x 10" 2 


E 


1.3 


-5.9 x 10~ 8 


-6.1 x 10~ 9 


F 


7.3 


1.1 x 10" 1 


8.2 x 10~ 3 


G 


7.1 


2.1 x 10~ 2 


-5.5 x 10~ 2 


J 


21 


-1.2 


-0.12 


K 


57 


-6.3 x 10~ 4 


9.1 x 10~ 4 
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Figure 13. Magnetic dipole moment dio/R^ (bottom) and mag- 
netic quadrupole moment d 2 \/R% l (top, middle) for model D, 
normalised to the initial value of dig/R^ = 5.29 X 10 _7 Soi as 
a function of time (in units of the Alfven time scale). All other 
components vanish due to symmetry. 



mum of 5.5 times the initial value. Likewise, the quadrupole 
(fei peaks during the transition phase before settling down 
to a constant value. The final field is highly axisymmetric, 
deviating from perfect symmetry by \d,2i\Rm/dio = 0.8 per 
cent. The asymptotic values of di m for models A-K are listed 
in Table El 



4.3 Boundary conditions 

We perform a simulation (model E) with the same initial 
configuration as model D (M a /M c = 1.0 and 6 = 3) but 
with inflow boundary conditions at R = R m . Fig. [T3] com- 
pares the density (left panel) and magnetic field (right panel) 
of models D and E. Again, inflow pins the magnetic field at 
the outer boundary, as opposed to outflow, which leaves the 
field free. The density distribution is almost unaffected. The 
magnetic field is mainly affected in the outermost region, 
where the plasma density is low. The overall time evolu- 
tion (a nonaxisymmetric transition phase which leads to a 
nearly axisymmetric equilibrium) remains as before, too. We 
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Figure 14. Meridional section of density contours (left) and mag- 
netic field lines (right) for model D (solid curve) and model E (dot- 
ted curve). The inflow boundary condition corresponds to line- 
tying at the outer boundary. Deviations between the two models 
occur in the outermost, low density regions. 



therefore conclude that the outer boundary condition can be 
chosen opportunistically. 

By contrast, the inner boundary condition contributes 
fundamentally to stability. The tension of the magnetic field, 
which is tied to the stellar surface, suppresses those modes 
which are driven by a pressure gradient perpendicular to 
the magnetic flux surfaces, such as the interchange and bal- 
looning mode. If line-tying is taken away, the latter modes 
disrupt the mountain in short order. If we rerun model A (for 
example) by applying a reflecting boundary condition at 
r = R», the mountain rapidly dissolves on a timescale ~ To ■ 
The same experiment for model D results in high velocities 
and steep field gradients, causing the numerical algorithm 
of zeus-mp to break down. 



4.4 Energetics 

iMouschoviasI (|l974l ) showed that an isothermal, gravitating, 
MHD system possesses a total energy W , which can be writ- 
ten as the sum of gravitational, kinetic, magnetic, and acous- 
tic contributions, defined by the following volume integrals, 
evaluated over the simulation volume: 



W s = / dV PV , 



(10) 
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Figure 15. The evolution of the total energy W and its com- 
ponents Wm, W g , Wk, and W a (top to bottom) for model D, all 
normalised to Wo = 2.2 X 10 36 erg, as a function of time (in units 
of the Alfven time scale). The total energy increases artifically, 
due to mass loss through the outer boundary (see text). 



Figure 16. The evolution of the total energy W and its com- 
ponents Wm, Wg, Wk, and Wa (top to bottom) for model D, all 
normalised to Wo = 2.2 X 10 36 erg, as a function of time (in units 
of the Alfven time scale). W, W g , and W a are now corrected for 
the mass loss through the outer boundary (cf. Fig. 115 1 - 



W k = ~ J dVpv 2 , (11) 
W m = J- / dVB 2 , (12) 

W*= [ dVplogp. (13) 

Here, v is the plasma velocity and p = c 2 p is the pressure. 

The evolution of (|l(J[l - (|13p for model D is shown in Fig. 
1151 The magnetic energy (second panel from top) steadily 
decreases to 20 per cent of its original value, as the axisym- 
metric equilibrium evolves to a lower energy, nonaxisym- 
metric state. The kinetic energy peaks at t — 1.2ta, during 
the transition phase when the magnetic reconfiguration oc- 
curs. However, the gravitational and acoustic contributions, 
which dominate W, increase with time. The reason for this 
becomes apparent if we track the total mass in the simu- 
lation volume. Approximately 3.7 per cent of the mass is 
lost through the outflow boundary at r = R m by t = 6ta- 
The mass loss is responsible for the increase of W a oc p 2 and 
Wg oc p, both of which are negative (W g because the plasma 
is gravitationally bound and W a since p < 1 in our units). 

Let us try to correct for the mass loss by multiplying 
W s , W fe , and W a by M(t = 0)/M(t), where M(t) is the mass 
in the simulation volume at time t. The result is presented 
in Fig. [16] W 9 and W a now decrease, and the total energy, 
W — W g + Wk + Wm + W a , decreases by just 2.5 per cent. 

The approximate correction above assumes p decreases 
uniformly, which is not strictly true. We therefore check our 
claim that mass loss is responsible by tracking the energy 
evolution of model E, which has the same initial configu- 
ration as model D, but an inflow outer boundary which 
blocks mass loss. From Fig. 1171 it is clear that the total en- 
ergy rises then falls, consistent with the observed dynamical 
evolution. The mountain oscillates until toroidal modes grow 
sufficiently to disrupt the initial configuration and force it 
into a nonaxisymmetric state. There is no spurious increase 
in W. We conclude that mass loss through the outer bound- 




Figure 17. The evolution of the total energy W and its com- 
ponents W m , W g , VKk, and W& (top to bottom) for model E, all 
normalised to Wo = 2.2 X 10 36 erg, as a function of time (in units 
of the Alfven time scale). 

ary is indeed responsible for the observed behaviour of W 
in model D in Fig. 1151 

4.5 Dependence on AI a 

Does the final, nonaxisymmetric configuration of the moun- 
tain become unstable once the accreted mass exceeds a crit- 
ical threshold? There are two ways that this can happen. 
First, the sequence of nonaxisymmetric GS equilibria passed 
through as M a increases can terminate above a critical value 
of M a \ i.e. there is a loss of equilibrium. PM07 observed 
this phenomenon in axisymmetric magnetic mountains with 
M a > 1CP 4 Mq, when the source term in the GS equation 
forces the flux function outside the range ^ %p ^ ip* per- 
mitted by the boundary condition at r = R*. Second, the 
nonaxisymmetric state reached in Fig. [7] (for example) may 
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Figure 18. The mass quadrupolc moments for models F (solid), 
D (dotted), and G (dashed), normalised to 1.33 X 10 25 g cm 2 , 
as a function of time. While all models show similar dynamical 
behaviour, the quadrupole moment of the final state increases 
with M a . 



Figure 19. Ellipticity e for models F (solid), D (dotted), and G 
(dashed) as a function of time. As expected, e increases with M a . 



3.5 



be metastable. That is, it may be a local energy minimum 
which can be reached from an axisymmetric starting point 
via the Parker instability but which the system can exit (in 
favor of some other, global energy minimum) if the system 
is kicked hard enough. One way to kick the system hard is 
to increase M a substantially. 

We are not really in a position to answer this ques- 
tion definitively, because the GS fails to converge to valid 
equilibria for M a ^> 10~ 4 M o , due to numerical difficulties 
(steep gradients, which would be smoothed in a more real- 
istic, non-ideal-MHD simulation). Nevertheless, we begin to 
address the issue by performing two simulations, models F 
and G, with the same parameters as model D but with lower 
and higher masses viz. M a /M c = 0.6 and M a /M c — 1.4 re- 
spectively. The mass quadrupole moments are plotted versus 
time in Fig. 1181 The solid and dashed curves are for models F 
and G respectively, with model D (dotted curve) overplotted 
for comparison. 

The dynamical behaviour of all three models is simi- 
lar: a violent transition phase which settles down to a non- 
axisymmetric state. However, the start of the transition 
phase, defined as the instant where Q33 is maximized, scales 
roughly oc 0.5M a /M c in units of ta- Physically, this means 
that the onset of the toroidal instability depends on M a . 
We can understand the trend in terms of the Parker insta- 
bility (section I4.1[l . whose growth rate scales as Tp oc w^ 1 . 
By measuring va at 9 = 7r/2 in models D,F, and G, we find 
va oc M a empirically and Tp oc M~ , consistent with the 
Parker scalings. 

The evolution of the ellipticity e oc Q33 for models D, 
F, and G is displayed in Fig. 1191 Of chief interest here is 
the ellipticity o f the final state. It increa ses along with M a , 
consistent with iMelatos fe Pavnd (|2005l ). A linear fit yields 
the following rule of thumb for our downscaled star (section 
[13}: 



Note, however, that the fit is valid in the range 0.6 ^ 
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Figure 20. Azimuthal magnetic field component \B^,\ for models 
F (solid), D (dotted), and G (dashed), in units of Bo, plotted as 
a function of time, in units of the Alfven time. 

M a /M c ^ 1.4. Numerical di fficulties prevent u s from e xtend- 
ing it to larger values of M„.. |Pavne fc Melatosl (|2006al ) found 
e/10" 10 = 7.8M a /M c (l + l.lMa/Mc)" 1 in the M a ~ M c 
regime for the axisymmetric equilibrium. Equation (|14|l 
yields values roughly 70 per cent lower than the latter for- 
mula. 

For completeness, we plot the magnitude of the toroidal 
field component \B$\ versus time in Fig. 1201 Interestingly, 
the peak value is achieved for the intermediate mass model, 
D, not for model G. However, B^ in the final state depends 
weakly on M a . We find B^.f = 3.4 x 10~ 7 B , B^q = 4.5 x 
10 -7 -Bo, and B^,g = 3.2 x 10 _7 i?o, where Bo is defined in 
section 12.31 These values are comparable to the magnitude 
of the polar magnetic field B p = 2.9 x 10 _7 i?o. 

4.6 Dependence on curvature 

As discussed in section 12.31 PM07 argued that reducing R, 
and A4* does not affect the equilibrium structure as long as 
ho remains constant, at least in the small- M a limit. To test 
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Figure 21. Mass quadrupole moments Qij for model J (a = 75, 
solid curve) and model K (a = 100, dashed curve), plotted as 
a function of time (in units of the Alfven time). The scale for 
model J (K) appears on the left (right) vertical axis. Although 
the transition phase is less distinct than in Fig. 1111 these models 
basically share the same dynamics as lower curvature runs. 



whether this also holds for the dynamical behaviour of the 
system, we perform two runs, models J and K, with a — 75 
and a = 100 respectively. 

Fig- [2Dpl°ts Qij versus time for these models. The tran- 
sition phase is more gradual than model D (Fig. Again, 
however, Q12 rises significantly , mark ing a deviation from 
axisymmetry. iMelatos fc Payne! (120051 ) found £ tx a analyt- 
ically in the small-M a regime, so we fit a parabola to the 
simulation data (for M a = M c ): 



10- 



1.82a 



(15) 



A realistic star has a = 1.9 x 10 4 (cf. section l!0|) . Ex- 
trapolating (|15[) , we find e = 6.6 x 10 -5 . (An ellipticity 
this large is close to the upper limit inferred from exist- 
ing gravitational-wave nondetections; see section[6]for more 
details.) However, it should be remembered that equation 
()15[) is an overestimate, because the computations in this 
paper neglect nonideal MHD effects. 



5 GLOBAL MHD OSCILLATIONS 

In this section, we explore the natural oscillation modes of a 
nonaxisymmetric magnetic mountain. We do this by loading 
the final state from models D, F, and G into zeus-mp and 
setting v = on the whole grid. This procedure introduces 
numerical perturbations that are sufficient to excite small 
linear oscillation modes, albeit an uncontrolled distribution 
thereof. We then compute the power spectrum 



P[S\t 



J2 s ^ 



Jti/N 



(16) 



by evaluating the discrete Fourier transform of the scalar 
function S(t) [e.g. B r (t)] at sample times ti. 

In order to explore the magnetic modes, we examine 
B r , Be, and B^. We choose one point on the grid where 
the amplitude of the oscillations is high, namely (r, 6, (f>) = 



Figure 22. Power spectrum of B r , Bg, and B^ (top to bottom) 
for model D (arbitrary units), plotted as a function of angular 
Fourier frequency (in units of Tq 1 ) . 



(50.03/io,0.26,2.06) and compute P[B r ], P[B e ], and P[B^}. 
The results are displayed in Fig. 1221 We can distinguish five 
different spectral peaks at 10 3 wt = 0.78,1.1,1.2,2.1,3.2, 
which are more or less distinct for the different components. 

For a magnetized gravitating slab in a plane-parallel 
geometry, one can distingui sh three different MHD modes 
jGoedbloed fc Poedtsl |2004| ): slow magnetosonic, Alfven, 
and fast magnetosonic. Each mode consists of a discrete set 
of eigenmodes and a continuous spectrum, which are clearly 
separated. Unfortunately, such clean separation cannot be 
expected for a highly inhomogenous plasma in spherical ge- 
ometry. Generally, different parts of the spectrum overlap 
or degenerate into a single point in a nontrivial way. We 
therefore restrict the discussion below to some qualitative 
remarks. 

The MHD spectrum contains genuine singularities, 
when the eigenfrequency coincides with the Alfven or slow 
magneto-sonic frequency at some location within the mag- 
netic mountain. In this case, the boundary value problem be- 
comes singular; the boundary conditions can be fulfilled for 
a continuous range of frequencies. The singular frequencies 
depend on the components of the wave vector perpendicular 
to the direction of inhomogenity. 

It is unclear whether the band uj < 0.002^ , which 
looks "filled" in Fig. 1221 belongs to the continuous part of 
the spectrum or else is an artifact of the nonzero line width 
from numerical damping (which can be estimated from the 
sample times U to be ~ 1.2 x 10~ 4 ojto). We do not ob- 
serve any singular behaviour in the field variables, but we 
note that singularities would be suppressed by the shock- 
capturing algorithm (i.e. the artificial viscosity) in zeus-mp. 
We conclude that the features in Fig. [22] are probably dis- 
crete lines. 

Let us compare these results to the spectrum of the 
axisymmetric model A (Fig. I23|l . We first note that the 
Alfven frequency cja = O-OlSr^ 1 (M a /M c ) -1 ^ 2 and acous- 
tic frequency uj s = 0.48r _1 found by PM07 are outside 
the range of this plot, which is set by the Nyquist fre- 
quency un = (4-ivNAt)~ (At — 50to for model A and 
At = IOtq for models D-G). Here, we are restricted to 
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A complete analytic determination of the discrete and 
continuous components of the MHD spectrum via a full lin- 
ear mode analysis will be attemped in a forthcoming paper. 
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Figure 23. Power spectra of B r and Bg (top to bottom) for 
model A (arbitrary units), plotted as a function of angular Fourier 
frequency (in units of rj~ ). 




Figure 24. Power spectra of B r , Bg, and B^ (top to bottom) for 
models F (solid curve) and G (dashed curve) in arbitrary units, 
plotted as a function of angular Fourier frequency in units of Tq . 
We overplot the spectrum of model G, stretched by a factor of 
1.4 in ui, as dotted curve, by way of comparison. 



low frequency oscillations which are generally associated 
with global magnetic modes. Most distinct is the peak at 
10 3 ojto = 0.5, which is not visible in Fig. [22] This long 
wavelength poloidal mode is suppressed in favor of toroidal 
modes in the three-dimensional configuration. However, the 
small peak at 10 3 ljto = 1.1 is present in both systems. This 
example illustrates vividly how relaxing the axisymmetric 
constraint leads to a different MHD spectrum. 

Fig. 1241 shows the power spectrum for models F (solid) 
and G (dashed). The most distinct peaks are again concen- 
trated in the low frequency region. We can roughly match 
the peaks of models G and F by stretching the former spec- 
trum by a factor of 1.4 in frequency. The higher M a equi- 
librium has a similar structure, but the Alfven timescale is 
lower because the plasma density is 85 per cent higher. 



6 DISCUSSION 

Magnetically confined mountains on accreting neutron stars 
screen the magnetic dipole moment of the star. Potentially, 
therefore, the process of polar magnetic burial can explain 
the observed reduction of /j with M a in neutron stars with 
an accretion history. However, before magnetic burial can 
be invoked as a viable explanation, the question of stability 
must be resolved. In this article, we concentrate on the im- 
portant aspect of three-dimensional stability, deferring resis- 
tive processes to future work (especially the issue of resistive 
g-mode^f|). 

We find that the axisymmetric configurations in PM04 
are susceptible to the three-dimensional magnetic buoyancy 
instability. The instability proceeds via the undular sub- 
mode, with growth rate oc A 1 / 2 , limited by the toroidal grid 
resolution. However, instead of breaking up and reverting 
to an isothermal atmosphere threaded by a dipolar mag- 
netic field, the magnetic field reconfigures (over a few Alfven 
times) and settles down into a new nonaxisymmetric equilib- 
rium which is still highly distorted. Just as the axisymmetric 
solutions in PM04 are the final saturated states of the non- 
linear evolution of the Parker instability in two dimensions, 
we find here the three-dimensional equivalent. This surpris- 
ing result is the main conclusion of the paper. It holds ir- 
respective of the outer boundary condition and curvature 
rescaling factor, but it depends critically on the line-tying 
boundary condition at the stellar surface. 

The final state is predominantly axisymmetric, with 

1.5 < IQ12/Q33I/IO" 3 < 3.2 for models D, F, and G 
(0.6 < Ma/M c < 1.4). The ellipticity for model G reaches 

1.6 x 10~ 10 in the downscaled star. 

The stability of magnetic mountains is important for 
the emission of gravitational waves from accreting millisec- 
ond p ulsars, as pointed out previously bv lMelatos fc Pavnel 
(|2005l ). Persistent X -ray pulsations from accreting bi- 
nary pulsars imply that the angle between the spin vec- 
tor O and the magnetic symmetry axis p, is n ot zero 
(IRomanova et al ]|2004l ; iKulkarni fc Romanova|[2005l ) . Hence 
a magnetic mountain constitutes a time-varying mass 
quadrupole which emits gravitational waves. Furthermore, 
the star precesses in general, emitting gravitational waves 
at the spin frequency and its first harmonic. The amplitude 
of the resulting signal (with curvature upscaled to a realistic 
neutron star at a distance d = 10 kpc using e oc a 2 ) is plot- 
ted in Fig. [25] for 10~ 9 ^ M a /M Q sC 10" 3 . The amplitude 
of the average signal that can be detected by the Laser In- 
terferometer Gravitational Wave Observatory (LIGO) from 
a periodic source with a false alarm rate of 1 per cent and a 
false dismissal rat e of 10 per cent over an integration time 
of Tp = 14 days (jjaranowski et all Il998l ; I Abbott. B. et all 
120041 ). is overplotted in Fig. 1251 This To can realistically be 
achieved computationally. 

At this point, it is important to acknowledge that M a = 
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Figure 25. Amplitude of the gravitational wave signal ho 
for M a /M Q = 10- 9 ,10- 8 ,10" 7 ,10~ 6 ,10~ 5 ,l(r 4 ,l(r 3 , for 
the axisymmctric (solid lines) and nonaxisymmetric equilibrium 
(dashed lines). The sensitivities of Initial and Advanced LIGO, as- 
suming 14 days coherent integration, are also plotted (upper and 
lower curves respectively) . The growth of the mountain is arrested 
for M a > 1.2 X 10~ 5 M(T) (light shaded region), due to Ohm ic dis- 
sipation lMelatos fe Pavndl2005l ; lvigelius fe Melatosll2008f) . while 
the right-hand edge is excluded at present because no accreting 
millisecond pulsars have been discovered with /„ > 0.7 kHz (dark 
shaded region). 



1.4M C w 1.7 x 10" 4 M Q is still well below M a ~ O.IMq, 
the mass req uired to spin up a neutron star to millisec- 
ond periods l|Burderi et alj [l999). At present, this high- 
mass regime is not accessible numerically; neither the GS 
solver nor zeus-mp can handle the steep magnetic gradi- 
ents involved. By the same token, Ohmic diffusion becomes 
important in this high-M a regime (|Melatos fc Payne] 120051 ; 
IVigelius fe MelatoslfeOOct ). smoothing the gradients and mit- 
igating the numerical challenge. We postpone studying real- 
istic values of M a to an accompanying paper, which will con- 
centrate on non-ideal MHD simulations. However, to make 
a rough estimate regarding detectability here, we assume 
that non-ideal effec ts stall the growth of th e mountain at 
M a « M c , following iMelatos fe PaTnel i|2005l ). This includes 
the region shaded light grey in Fig. [55] Furthermore, no ac- 
creting millisecond pulsars have been discovered spinning 
faster than f t > 720 Hz, possibly due to braking by gr avita- 
tional waves (|Bildstenlll998l ; IChakrabartv et alj|200ot ). The 
region with 2/* > 1400 Hz is shaded dark grey in Fig. 1251 

Even with those exclusions, Fig. [25] demonstrates that 
there is a fair prospect of detecting gravitational waves 
from accreting X-ray millisecond pulsars in the near future, 
for accreted masses as low as M a ~ 10~ 4 Mq. Recent di- 
rected searches for gravitational waves from the nearby X- 
ray source Sco-Xl foun d no signal at the level ho > 10~ 22 
|Abbott. B. et al.ll2007l ) , thereby setting an upper bound on 
the ellipticity of e = 3.6 x 10" 3 . 
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APPENDIX A: DEFINING THE GRID AND 
BOUNDARY CONDITIONS IN ZEUS-MP 

In this appendix, we briefly outline the key variables and 
settings in zeus-mp, to aid the reader in reproducing 
our numerical results. Our grid consists of ggenl:nbxl, 
ggen2:nbxl, and ggen3:nbxl blocks in the r, 9, and <j> 
direction, respectively. The integration volume is defined 
by FU/ho < r < Rm, 0^9^ tt/2, ^ (j> < 2tt. 
The radial coordinate in the GS code, x = (r — R»)/ho 
(0 ^ x ^ X), is stretched logarithmically according to 
xi — log(:r + e~ La: ) + L x , where L x controls the zooming 
(PM04). This grid is implemented by setting the zeus-mp 
parameters ggenl : xlmin — R*/ho and ggenl : xlmax = 
ggenl : xlmin + X. Stretching is achieved via the param- 
eter ggenl :xlrat, which sets the radial length ratio of two 
neighbouring zones. In order to get consistent radial grid po- 
sitions in the GS code and zeus-mp, we set ggenl : xlrat = 
(Xe L * +1) (G *- 1)-1 . 

Boundary conditions are enforced in zeus-mp via ghost 
cells, which frame the active grid cells. Several predefined 
prescriptions are supplied to implement a variety of standard 
boundary conditions. In the <j> direction, we choose periodic 
boundary conditions [ikb.niks(l)=4 and okb.noks(l)=4]. 
The 9 = 7r/2 surface is reflecting, with normal magnetic 
field [ojb.nojs(l)= 5], which translates to vi = B|| = 0. 
The line 9 = is also reflecting [ijb.nijs(l)= -1] with 
tangential magnetic field (v_i_ = Bi = 0). Additionally, 
the toroidal component is reversed at the boundary, i.e. 
B^ = —B^, where B^ and B^ are the field components 
for 8 < and 9 > 0, respectively. The outer surface r = R m 
is usually an outflow [oib nois(l)= 2] boundary, i.e. zero 
gradient. The stellar surface is impenetrable, so the inner 
r = R, boundary is inflow [iib.niis(l)= 3]. This enables 
us to impose line-tying at r = R, by fixing the density and 
magnetic field there. We also use an isothermal equation of 
state (XIS0=.true.). 



APPENDIX B: MASS MULTIPOLE MOMENTS 

We work out the mass quadrupole moment in Cartesian co- 
ordina tes from the cod e output in spherical coordinates. Fol- 
lowing |jackson| (| 199Sf ) . we define the spherical mass multi- 
pole moments according to 



d 3 x'ir m (0',</>V P (x'), 



(Bl) 



where Y; m denotes the usual orthonormal set of spherical 
harmonics. 

The spherical quadrupole moments are related to the 
traceless, Cartesian quadrupole moment tensor, 



Qi 



d 3 x' {"ix'iXj — r' 2 Sij)p(x.'), 



by Qn 



6(27r/15) 1/2 Re(g : 



(B2) 

2(47r/5) 1/2 g 2 o, 

1/2t 



Q12 = -6(27r/15) 1/2 Im(g 22 ), Q13 = -3(87r/15) 1/2 Re(g 2 i), 
Q 22 = -6(27r/15) 1/2 Re( g22 ) - (47r/5) 1/2 g 20 , and 
Q23 = 3(87r/15) 1/2 Im( 92 i). 

In the axisymmetric case (when the star and the moun- 
tain form a prolate spheroid), we have p = p(r, 9) and the <j> 
integrals in (|B1|I vanish. Q is then diagonal with components 
Qn = Qyy = — Qzz/2, with respect to the body coordinate 
system, and we can introduce the ellipticitjQ e, where we 
assume that the z axis is the symmetry axis: 



3(/if 



la 



(B3) 



with la = 2M„i? 2 /5 is the moment of inertia along the 
rotation axis for a biaxial ellipsoid with mass M* and minor 
axis R*. We can compute the ellipticity directly from the 
code output through 



2h 



d9 1 



idr pr 4 sin 9(3 cos 2 9 — 1). 



(B4) 



APPENDIX C: MAGNETIC MULTIPOLE 
MOMENTS 

In a source- free region J = 0, a magnetic field 
B is determined solely b y its radial component B T 
l|Bouwkamp fc Casimir|[l954h . which, from Maxwell's equa- 
tions, satisfies the Laplace equation 



V 2 B r = 0. 



(CI) 



One can therefore define the magnetic multipoles as the ex- 
pansion coefficients in the general solution of the boundary 
value problem (|C1[1 . viz. 



;=o m=—l 

with 



kr n r 



(C2) 



4 Note that the definition of the ellipticity varies in the 
literature. The ellipticity define d here is consistent with 
iBonazzola &: Gourgoulhorj l ll996l) an d ISh apiro fc Teuko lskvl 
(Il983ll and is related t o the e llipticity in lMelatos fc Payne! (1200511 
and [ jaranowski et al. I GUI) by |e| = 3eMP- lAbbott. B. et al,l 
l|2007F used a different ellipticity defined for a triaxial rotator, 
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d lm = r l+1 J AO, Y t * m r ■ B. (C3) 

Note that dio is related to the magnetic moment fi of 
a dipole field B(r) = ^r- _3 (2cos^e r + sin 8ee) by dio = 
4( 7 r/3) V V- 

In the case of north-south symmetry, we find dio = 2dio 
and e?2i = 2<i2i, where a hat denotes the moment evaluated 
on the hemisphere. All other coefficients vanish. 



